import libnum


def getFullP(low_p, n):
    R. < x > = PolynomialRing(Zmod(n), implementation='NTL')
    p = x * 2 ^ bit + low_p
    root = (p - n).monic().small_roots(X=2 ^ 128, beta=0.4)
    if root:
        return p(root[0])
    return None


def phase4(low_d, n, c):
    maybe_p = []
    for k in range(1, 4):
        p = var('p')
        p0 = solve_mod([3 * p * low_d == p + k * (n * p - p ^ 2 - n + p)], 2 ^ bit)
        maybe_p += [int(x[0]) for x in p0]
    # print(maybe_p)

    for x in maybe_p:
        P = getFullP(x, n)
        if P: break

    P = int(P)
    Q = n // P

    assert P * Q == n
    print(P)
    print(Q)

    d = inverse_mod(3, (P - 1) * (Q - 1))
    print(d)
    print(libnum.n2s(int(power_mod(c, d, n))))


bit = 486
n = 99233273001596380809501393613886417854988989363311895592445631124571940638105064581031336422148895730117448547059137492090852484233539441265373247940823283633017582469362503632785297924194187912199716752955609363457416190782142095008241313065484612071705711161086869849047664563205898255359661312582650481473
e = 3
c = 175676150266403937224898870626869248307097859453341599800113943191154294552011908698393750389195590199207971365632903719917006078351629939912360175671032635640354766675409868021903917260597989036476083685690139071290022752606720020238530580507331902326179201548484337939338062870542095137303322195300197
d1 = 66155515334397587206334262409257611903325992908874597061630420749714627092070043054020890948099263820078299031372758328060568322822359627510248831960548842374360111368587702677016263988651746874367056185352095585037936034196599496577057035033718066207964292781685419156037246789689219471468151618592221802699

phase4(d1, n, c)